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ABSTRACT 

Prototype  devices  that  use  pulse-power  techniques  to  generate  an  intense  acousti¬ 
cal  field  in  water  have  fostered  a  renewed  interest  in  applying  finite-amplitude  sound 
to  the  mine-neutralization  problem,  particularly  in  a  littoral  region.  A  simple  intu¬ 
itive  description  of  acoustical  mine  neutralization  includes  three  basic  processes:  (1) 
generation  of  the  acoustical  field  at  the  source;  (2)  nonlinear  wave  propagation;  and 
(3)  neutralization  mechanisms  at  the  target.  This  document  focuses  on  the  second 
issue,  the  propagation  of  an  intense  acoustical  field  from  a  source  to  a  target.  The 
research  discussed  here  provides  a  theoretical  foundation  for  a  modeling  effort,  de¬ 
scribes  several  case  studies,  and  gives  empirical  rules  for  establishing  the  material 
parameters  of  fresh  and  sea  water  required  by  the  theory.  Several  key  issues  are 
presented,  including  the  breakdown  of  using  linear  superposition  of  the  results  for 
discrete  sources  in  an  array,  phasing  an  array  of  discrete  sources  for  beam  steering, 
and  peak  positive  and  negative  pressures. 


TR  9803  iii 


UNIVERSITY  OF  WASHINGTON  •  APPLIED  PHYSICS  LABORATORY 


CONTENTS 

1.  INTRODUCTION  .  1 

2.  CASE  STUDIES  .  5 

2.1  Boston  University  Experiments  .  5 

2.2  Fresh  Water  Lake  .  11 

3.  THEORETICAL  DEVELOPMENT  .  15 

3.1  Augmented  Khokhlov-Zabolotskaya-Kuznetsov  Equation  .  15 

3.2  Axisymmetric  Circular  Source  .  15 

3.3  Rectangular  Source  .  18 

4.  NUMERICAL  SOLUTION  .  19 

4.1  Diffraction  Algorithm  .  20 

4.2  Thermoviscous  Absorption  Algorithm  .  21 

4.3  Quadratic  Nonlinearity  Algorithm  .  22 

4.4  Relaxation  Algorithm  .  23 

APPENDIX,  Empirical  Material  Parameters  for  Water  .  24 

REFERENCES  .  28 


iv  TR  9803 


UNIVERSITY  OF  WASHINGTON  •  APPLIED  PHYSICS  LABORATORY 


LIST  OF  FIGURES 

Figure  1.  Retarded  time  traces  along  the  acoustic  axis  at  25  kPa  .  7 

Figure  2.  Retarded  time  traces  along  the  acoustic  axis  at  125  kPa  .  7 

Figure  3.  Retarded  time  traces  along  the  acoustic  axis  at  225  kPa  .  8 

Figure  4.  Retarded  time  traces  along  the  acoustic  axis  at  325  kPa  .  8 

Figure  5.  Center  two  cycles  from  Figures  l(e)-4(e)  .  9 

Figure  6.  Time-averaged  energy  flux  along  acoustics  axis  .  10 

Figure  7.  Retarded  time  traces  produced  by  a  phased  array  in  a  fresh 

water  lake  .  13 

Figure  8.  Retarded  time  traces  for  a  phased  array  and  focused  piston 

transducer  .  13 

Figure  9.  Time-averaged  energy  flux  for  an  array  in  a  fresh  water  lake  .  14 

Figure  10.  Graphical  depiction  of  the  nonlinear  distortion  algorithm  .  22 


TR  9803  v 


UNIVERSITY  OF  WASHINGTON  •  APPLIED  PHYSICS  LABORATORY 


1.  INTRODUCTION 


Before  intense  acoustic  fields  can  be  employed  as  a  means  of  neutralizing  mines 
remotely,  several  critical  issues  regarding  nonlinear  wave  propagation  need  to  be  ad¬ 
dressed  and  clarified.  When  a  finite-amplitude  wave  propagates  in  a  fluid,  four  basic 
physical  phenomena  affect  the  wave  propagation.  These  phenomena  are  nonlinearity, 
thermoviscous  absorption,  relaxation,  and  diffraction.  A  brief  description  of  each  is 
given  below.  Other  critical  issues  for  the  neutralization  of  mines  with  finite-amplitude 
waves  include  environmental  conditions,  beam  forming  (phasing  an  array  of  discrete 
sources),  expected  peak  pressure  amplitude  (acoustic  saturation),  and  signal  dura¬ 
tion.  To  date,  no  comprehensive  study  has  considered  these  issues  as  they  relate 
to  nonlinear  wave  propagation  in  an  ocean  environment  and  a  mine  countermeasure 
(MCM)  scenario. 

Nonlinearity  has  important  implications  for  finite-amplitude  waves.  Nonlinear¬ 
ity  causes  steepening  of  the  pulse,  which  ultimately  leads  to  shock  formation.  The 
propagation  velocity  of  the  finite- amplitude  wave  depends  on  the  local  sound  pres¬ 
sure.  A  higher  positive  pressure  produces  a  higher  speed.  Thus,  a  wave  crest  moves 
faster  than  a  trough,  causing  the  wave  to  steepen.  As  a  consequence,  energy  from 
lower-frequency  components  is  repartitioned  to  higher-frequency  components.  Thus 
nonlinearity  is  not  a  direct  mechanism  of  energy  loss;  however,  the  pumping  of  en¬ 
ergy  to  higher  frequencies  leads  to  energy  loss  through  thermoviscous  absorption  and 
relaxation  processes. 

Energy  loss  through  thermoviscous  absorption  results  in  a  heating  of  the  fluid. 
A  finite-amplitude  wave  imparts  some  momentum  to  a  “parcel”  of  fluid.  As  adja¬ 
cent  parcels  of  fluid  “rub”  against  one  another,  energy  is  converted  to  heat  through 
frictional  forces.  This  mechanism  is  directly  related  to  the  thermal  and  viscose  prop¬ 
erties  of  the  fluid  through  the  attenuation  coefficient.  This  coefficient  is  frequency 
dependent  such  that  higher-frequency  components  of  a  pulse  will  be  more  strongly 
attenuated  than  lower-frequency  components.  For  fresh  and  sea  water,  the  thermo¬ 
viscous  attenuation  coefficient  goes  as  /2,  where  /  denotes  frequency. 

Although  relaxation  is  an  energy  loss  mechanism,  it  differs  from  thermoviscous 
absorption  because  it  depends  on  the  presence  of  impurities.  In  the  ocean,  the  two 
primary  sources  of  relaxation  are  magnesium  sulfate  [1],  MgS(>4,  and  boric  acid  [2], 
H3BO3,  which  dissipate  energy  through  chemical  activity  (e.g.,  chemical  dissociation 
or  increased  vibrational  or  rotational  energy  of  a  molecule).  In  addition,  relaxation 
processes  give  rise  to  dispersion  (frequency-dependent  phase  velocity)  which  leads  to 
temporal  spreading  of  the  pulse.  Higher-frequency  components  of  a  pulse  typically 
propagate  with  a  faster  phase  velocity,  stretching  the  pulse. 
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Diffraction  is  associated  with  the  finite  aperture  of  a  directive  source.  Its  main 
contribution  to  the  wave  propagation  problem  is  a  spatial  spreading  of  the  energy. 
Although  geometric  focusing,  beam  forming,  and  amplitude  shading  can  reduce  the 
spreading  for  a  bounded  beam  source,  diffraction  eventually  sets  in  beyond  the  focal 
zone. 

The  environment  has  severe  implications  for  finite-amplitude-wave  neutralization 
of  mines.  Bubbles  are  present  in  sufficient  quantity  in  shallow  and  very  shallow  water 
such  that  even  linear  acoustics  is  altered  [3,4].  These  bubbles  may  be  entrained 
from  breaking  waves  in  the  surf  zone  or  advected  from  deeper  water  by  currents. 
Bubbles  can  be  stabilized  by  a  thin  coating  of  surfactant,  and  the  lifetime  of  such  a 
bubble  can  greatly  exceed  the  lifetime  of  a  similar,  clean  bubble  in  fresh  water.  With 
respect  to  nonlinear  wave  propagation,  bubbles  cause  a  significant  increase  in  f3,  the 
coefficient  of  nonlinearity.  For  fresh  and  sea  water  without  bubbles,  /?  is  nominally 
3.5  and  3.6,  respectfully.  The  addition  of  even  a  small  number  density  and  bubble  size 
distribution  can  increase  /?  by  several  orders  of  magnitude,  with  values  reported  as 
large  as  31,000  [5].  With  this  increase  in  nonlinearity,  one  might  assume  that  finite- 
amplitude  wave  propagation  will  become  “easier”;  however,  the  bubbles  affect  the 
dispersion  and  attenuation  of  sound  in  water  as  well.  Thus,  the  presence  of  bubbles 
complicates  the  analysis  such  that  they  cannot  be  ignored. 

Bubbles  in  the  presence  of  an  intense  acoustic  field  provide  sites  for  cavitation. 
For  a  transient-type  shock  pulse,  a  region  of  negative  pressure  (relative  to  the  ambient 
hydrostatic  value)  follows  the  large  compressional  shock.  During  this  negative  pres¬ 
sure,  any  bubbles  present  will  experience  rapid  growth  from  their  initial  equilibrium 
size.  After  the  negative  pressure  passes,  the  bubbles  will  collapse,  with  subsequent 
acoustic  radiation.  This  radiation  may  be  beneficial  or  detrimental  to  MCM,  but 
proper  analyses  is  beyond  the  scope  of  this  document. 

Nonlinear  wave  propagation  in  shallow  and  very  shallow  water  necessarily  implies 
that  one  must  consider  the  air-sea  and  sea-sediment  boundaries.  For  a  transient 
excitation,  ray  tracing  provides  sufficient  information  to  determine  if  multipaths  will 
(or  will  not)  affect  the  direct  arrival  at  a  target.  For  example,  if  a  source  and  target  are 
5  m  below  the  air-sea  interface  and  separated  by  100  m,  then  the  first  surface  bounce 
will  have  a  0.5  m  longer  pathlength  than  the  direct  ray  in  a  homogeneous  ocean.  This 
implies  that  the  surface  bounce  will  arrive  at  the  target  approximately  0.3  ms  after 
the  direct  pulse.  If  the  pulse  duration  exceeds  0.3  ms,  then  the  direct  arrival  and  the 
multipath  arrival  will  interfere.  Furthermore,  the  reflection  of  a  pressure  pulse  from 
the  air-sea  interface  will  invert  the  phase  of  the  pulse.  Water  cannot  sustain  a  large 
negative  pressure  (relative  to  the  hydrostatic  pressure)  without  cavitation.  Hence, 
the  negative  portion  of  a  pressure  pulse  may  cause  a  local  increase  in  the  bubble 
number  density,  leading  to  additional  dissipation  near  the  interface. 
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When  linear  acoustics  applies,  beam  forming  is  a  well-known  method  for  focusing 
and  directing  acoustic  energy.  With  properly  designed  sources,  the  air-sea  and  sea- 
sediment  boundaries  are  no  longer  an  issue.  Unfortunately,  when  the  source  levels  of 
individual  elements  of  an  array  are  sufficient  for  nonlinearity  to  become  an  issue,  the 
degree  to  which  beam  forming  can  be  achieved  is  unknown.  Only  a  limited  investi¬ 
gation  has  been  conducted  of  the  nonlinear  behavior  of  two  discrete  point  sources  in 
the  quasilinear  approximation  [6,  7,  8].  A  generalization  of  that  analysis  to  an  array 
of  point  sources  appears  possible;  however,  an  array  of  discrete  spark-gap  sources, 
though  point-like,  will  produce  source  levels  that  invalidate  some  assumptions  re¬ 
quired  by  the  quasilinear  approximation.  Applying  ordinary  linear  acoustics  with 
absorption,  relaxation,  and  diffraction  to  a  finite-amplitude  wave  propagation  prob¬ 
lem  will  lead  to  an  erroneous  result.  Proper  inclusion  of  nonlinearity  is  essential  to 
correctly  model  nonlinear  wave  propagation  from  a  distribution  of  discrete,  intensive, 
sound  sources.  This  claim  will  be  substantiated  in  the  following  discussion. 

At  the  higher  source  levels  expected  of  pulse-power  discrete  arrays,  beam  forming 
may  become  more  difficult.  Consider  a  two-source  array  where  the  field  from  source  1 
is  time  delayed  relative  to  that  from  source  2,  and  the  fields  from  1  and  2  are  sufficient 
for  shock  formation.  The  Rankine-Hugoniot  conditions  for  source  1  suggest  that  its 
shock  front  will  propagate  with  a  higher  velocity  than  the  shock  front  from  source  2. 
This  implies  that  the  shock  front  from  source  1  may  overtake  the  shock  front  from 
source  2  prior  to  their  arrival  at  a  target  unless  the  initial  delay  properly  accounts 
for  the  higher  propagation  speed  of  the  shock  front  from  source  1.  The  subsequent 
interaction  of  the  fields  from  sources  1  and  2  may  not  produce  the  desired  effect  at 
the  target. 

The  peak  pressure  achievable  from  conventional  source  technology  appears  to  be 
limited  to  less  than  200  MPa  except  for  explosive  sources  (higher  peak  pressures 
may  be  achieved  from  focusing).  This  suggests  nonlinear  acoustics,  where  cubic  and 
higher-order  terms  in  the  standard  acoustic  field  variables  are  neglected,  can  be  used 
to  analyze  the  field  produced  by  these  sources.  Nonlinear  acoustics  applies  when  the 
acoustic  Mach  number  satisfies  e  =  po/poCo  ^  1>  where  p0,  Po>  and  c0  are  the  peak 
positive  pressure,  the  ambient  density  of  the  fluid,  and  the  small-signal  sound  speed, 
respectively.  The  definition  of  the  Mach  number  given  here  assumes  that  the  energy 
is  propagating  primarily  in  one  direction  and  a  progressive  plane  wave  approximation 
of  p0  =  poCqUq  holds,  where  u0  is  a  typical  acoustic  particle  velocity.  Hence,  a  peak 
pressure  of  200  MPa  gives  e  =  po/pocl  «  0.1,  and  the  theoretical  treatment  discussed 
here  can  be  used.  However,  if  the  peak  pressure  is  expected  to  exceed  200  MPa,  then 
nonlinear  acoustics  may  not  yield  reliable  results.  The  only  recourse  is  a  numerical 
solution  based  on  the  fundamental  equations  of  fluid  dynamics  (conservation  of  mass, 
momentum,  and  energy,  and  an  appropriate  equation  of  state  for  the  fluid). 
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Initial  computations  for  the  nonlinear  wave  propagation  can  be  based  on  either 
the  Khokhlov-Zabolotskaya-Kuznetsov  (KZK)  equation  or  the  nonlinear  progressive 
wave  equation  (NPE).  Both  the  KZK  and  NPE  are  parabolic  wave  equations  which  are 
valid  within  a  small  angle  along  the  axis  of  propagation.  Methods  for  the  numerical 
computation  of  the  KZK  or  NPE  model  are  easily  found  in  the  recent  acoustics 
literature  [9,  10,  11,  12].  An  augmented  KZK  equation  explicitly  includes  quadratic 
nonlinearity,  thermoviscous  absorption,  relaxation,  and  diffraction,  which  permits 
independent  investigation  of  each  physical  mechanism.  This  document  concentrates 
on  theoretical  results  produced  by  a  numerical  implementation  of  an  augmented  KZK 
equation  [13]. 

Acoustic  saturation  is  normally  associated  with  continuous  wave  (CW)  radiation 
from  a  source  at  a  fundamental  frequency  with  sufficient  intensity  to  cause  nonlinear 
processes  to  be  important  [14,  15].  However,  acoustic  saturation  may  affect  a  field 
from  a  source  with  a  finite,  but  long,  pulse  with  a  fundamental  frequency  (e.g.,  a 
tone  burst  from  a  magnetostrictive  piston  source).  A  simple  explanation  of  acoustic 
saturation  is  that  energy  in  the  fundamental  frequency  is  pumped  into  higher  harmon¬ 
ics.  That  is,  once  the  fundamental  frequency  component  has  reached  saturation,  any 
further  attempt  to  increase  the  energy  in  the  fundamental  frequency  dumps  energy 
into  the  higher  harmonics.  Thus,  the  pulse  undergoes  nonlinear  distortion,  leading 
to  shock  formation.  Acoustic  saturation  may  not  be  an  important  effect  for  the  tran¬ 
sient  fields  expected  from  pulse-power  systems  because  quadratic  nonlinearity  is  a 
cumulative  effect. 
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2.  CASE  STUDIES 

2.1  Boston  University  Experiments 

The  first  case  study  considers  conditions  that  are  well  suited  to  test  the  robustness 
of  the  theory  and  numerical  algorithms  discussed  in  greater  detail  in  Sections  3  and  4. 
Edson  and  Roy  at  Boston  University  (BU)  conducted  several  experiments  involving 
nonlinear  wave  propagation  from  an  eight-element  annular  array  [16].  Equipment 
failure  in  the  generation  of  eight  independent  source  waves  prevented  the  use  of  the 
outermost  annular  element  of  the  array,  and  thus  the  simulations  discussed  here  are 
limited  to  a  seven-element  annular  array. 

The  BU  experiments  were  conducted  in  filtered,  deionized  water.  The  dissolved 
gas  content  of  the  water  was  reduced  to  85%  of  saturation  to  inhibit  cavitation.  The 
ambient  pressure  was  approximately  1  atm,  and  the  salinity  of  the  water  throughout 
the  experiments  was  0  ppt  (i.e.,  fresh  water).  Attempts  to  maintain  a  constant  am¬ 
bient  temperature  were  not  implemented.  Each  experiment  started  with  an  ambient 
temperature  dictated  by  the  room  temperature  of  the  laboratory  (nominally,  17°C). 
Edson  and  Roy  reported  an  upward  drift  of  the  temperature  of  no  more  than  5°C 
over  the  duration  of  an  experiment.  They  attributed  the  drift  to  the  proximity  of  the 
water  tank  to  the  seven  power  amplifiers  used  to  drive  the  individual  elements  of  the 
array.  For  the  computations  included  here,  the  ambient  pressure,  temperature,  and 
salinity  are  set  to  1  atm,  17°C,  and  0  ppt.  All  required  material  parameters  for  the 
water  are  then  computed  via  the  empirical  rules  in  the  appendix. 

The  array  consisted  of  a  center  circular  element  and  six  concentric  rings.  Each 
element  was  constructed  from  a  1-3  piezocomposite,  the  center  frequency  of  the  array 
was  480  kHz,  and  the  bandwidth  was  210  kHz.  The  surface  area  of  an  element  was 
18.85  cm2,  and  the  radius  of  the  innermost  element  was  2.45  cm.  In  the  progressive 
wave  approximation  for  a  given  source  amplitude,  the  energy  flux  from  any  one  ele¬ 
ment  is  equivalent  to  that  from  any  of  the  others.  A  0.03-cm  gap  separated  adjacent 
elements,  which  were  coplanar.  The  elements  (and  auxiliary  electronics)  were  inde¬ 
pendent  and  isolated  to  reduce  cross-talk  and  permit  phasing  of  the  array  to  adjust 
the  depth  of  focus. 

The  source  waveform  reported  by  Edson  and  Roy  was  a  10-cycle  sine  wave  tone 
buist  with  temporal  shading.  No  amplitude  shading  across  the  array  was  used  in  the 
experiments.  The  pressure  in  the  source  plane  for  the  nth  element  was 

pn(R,  ct  =  0,t)=po  exp{— [(r  +  xpn)/Nn]2m}  sin (r  +  </>„),  (1) 

where  ipn  is  a  phase  relative  to  that  of  the  center  element  and  is  adjusted  to  focus 
the  array  at  24  cm  from  the  face  of  the  array.  Equation  (1)  is  expressed  in  a  form 
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such  that  R  and  a  are  dimensionless  radial  and  z  coordinates  and  r  is  a  dimensionless 
retarded  time.  The  coefficient  m  controls  the  rise  time  of  the  temporal  envelope  and 
N  is  the  number  of  cycles,  which  governs  the  total  length  of  the  pulse.  When  m  =  1, 
the  envelope  defined  by  the  exponential  factor  in  (1)  becomes  the  familiar  Gaussian 
envelope.  A  sequence  of  experiments  was  performed  in  which  p0  was  set  to  20,  44, 
200,  300,  400,  and  465  kPa.  For  the  computations  reported  here,  the  pulse  contains 
five’ cycles  and  p0  has  been  set  to  25,  125,  225,  and  325  kPa.  These  choices  were 
dictated  by  limitations  on  the  available  computational  facilities  (see  Section  4  for 
further  comments). 

In  the  computations,  the  phase  of  each  element  must  be  set  to  an  appropriate 
value  relative  to  the  phase  of  the  center  element  to  cause  the  desired  focusing.  This 
is  accomplished  by  setting  </>, 0  =  0  and  tpn  =  GR^(n  >  0),  where  G  is  the  focusing 
gain  (proportional  to  the  ratio  of  the  outer  radius  of  the  array  to  the  focal  length) 
and  R  is  the  mean  radius  of  an  element.  The  choice  of  R  is  somewhat  arbitrary;  it 
could  have  been  set,  for  example,  to  either  the  inner  or  outer  radius  of  each  ring. 

Figures  1-4  correspond  to  retarded  time  traces  at  various  locations,  a,  along 
the  acoustic  axis  of  the  array  for  source  amplitudes  of  25,  125,  225,  and  325  kPa, 
respectively.  These  locations  are  (a)  0.2,  (b)  0.4,  (c)  0.6,  (d)  0.8,  (e)  1.0,  and  (f)  1.2 
where  the  normalization  is  the  focal  length  of  24  cm.  Although  not  apparent  in  all  the 
figures,  panels  (a)-(d)  and  (f)  contain  two  curves,  a  solid  curve  showing  the  results 
when  all  array  elements  are  driven  simultaneously  with  the  proper  phase,  and  a  long 
dashed  curve  showing  a  linear  superposition  of  the  possibly  nonlinear  acoustic  field 
when  each  element  is  driven  independently  (but  with  the  proper  phase).  Panel  (e) 
in  each  figure  is  similar  to  the  other  panels  except  a  third  curve  is  included:  a  short 
dashed  curve  that  shows  the  time  trace  produced  by  a  single  continuous  element 
that  is  spherically  focused  at  24  cm  (i.e.,  a  focused  circular  piston  transducer).  The 
significance  of  this  curve  is  that  it  is  the  upper  limit  on  the  expected  field  when  the 
number  of  array  elements  becomes  large.  Finally,  as  a  side  note,  the  shape  of  the 
pulse  in  Figure  1(d)  is  essentially  the  shape  of  the  pulse  in  the  initial  source  plane  as 
given  by  (1). 

Clearly,  Figure  1  demonstrates  that  at  a  low  source  amplitude  linear  superposition 
may  be  applied  to  the  acoustic  fields  produced  by  an  array  of  discrete  sources.  This 
permits  a  substantial  reduction  in  the  computational  complexity  and  allows  more 
general  array  geometries.  As  the  source  amplitude  is  increased,  however,  the  effects  of 
nonlinearity  become  important,  particularly  in  the  vicinity  of  the  focus.  To  illustrate 
this  point,  the  center  two  cycles  of  Figures  1(e),  2(e),  3(e),  and  4(e)  are  reproduced  in 
Figure  5  without  the  additional  third  curve.  To  reiterate,  the  solid  curve  represents 
the  acoustic  field  when  all  elements  are  operated  simultaneously,  and  the  dashed 
curve  depicts  the  linear  superposition  of  the  fields  from  the  individual  elements  when 
nonlinear  wave  propagation  is  permitted  for  those  fields. 
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Figure  2.  Retarded  time  traces  along  the  acoustic  axis  of  a  seven-element  annular 
array  at  the  following  a  locations:  (a)  0.2,  (b)  0.4,  (c)  0.6,  (d)  0.8,  (e)  1.0,  and  (f) 
1.2.  The  pressure  amplitude  in  the  source  plane  is  po  =  125  kPa. 
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Figure  3-  Retarded  time  traces  along  the  acoustic  axis  of  a  seven-element  annular 
array  at  the  following  a  locations:  (a)  0.2,  (b)  0.4,  (c)  0.6,  (d)  0.8,  (e)  1.0,  and  (f) 
1.2.  The  pressure  amplitude  in  the  source  plane  is  p0  =  225  kPa. 


Figure  4.  Retarded  time  traces  along  the  acoustic  axis  of  a  seven-element  annular 
array  at  the  following  a  locations:  (a)  0.2,  (b)  0.4,  (c)  0.6,  (d)  0.8,  (e)  1.0,  and  (f) 
1.2.  The  pressure  amplitude  in  the  source  plane  is  po  =  325  kPa. 
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Figure  5.  The  center  two  cycles  from  (a)  Figure  1(e),  (b)  Figure  2(e),  (c)  Figure 
3(e),  and  (d)  Figure  4(e).  The  solid  curve  is  the  results  when  driving  the  array 
elements  simultaneously  (with  the  proper  phase).  The  dashed  curve  corresponds  to 
a  linear  superposition  of  the  acoustic  fields  from  the  individual  elements. 


Inspection  of  Figure  5  shows  that  three  important  features  occur  as  the  source 
amplitude  is  increased.  First,  the  expected  peak  positive  pressure,  P+,  is  underpre¬ 
dicted  by  linear  superposition  of  the  independent  nonlinear  fields.  At  the  highest 
source  level,  P+  from  a  properly  phased  array  is  approximately  65%  greater  than  the 
simple  linear  superposition  predicts.  The  second  feature  involves  the  peak  negative 
pressure,  P_.  Linear  superposition  overpredicts  P-  although  the  severity  of  the  error 
is  only  about  20%.  The  importance  of  P+  and  P_  with  respect  to  mine  neutralization 
has  not  been  established.  The  third  and  final  feature  is  the  shape  of  the  two  wave¬ 
forms  in  each  panel  of  Figure  5.  It  is  well  known  that  nonlinear  wave  propagation 
shifts  energy  from  the  fundamental  frequency,  500  kHz  for  this  case  study,  to  higher 
harmonics.  The  degree  of  waveform  distortion  is  directly  related  to  the  content  of 
the  higher  harmonics  produced  such  that  a  more  distorted  pulse  contains  more  har¬ 
monics.  When  the  acoustic  fields  are  propagated  independently  and  then  linearly 
summed,  the  shifting  of  energy  to  higher  harmonics  is  artificially  quenched. 

Although  the  source  amplitudes  used  to  produce  Figure  1-5  differ  from  those 
reported  by  Edson  and  Roy,  these  computations  support  their  conclusions  because 
their  source  levels  of  20,  100,  200,  and  300  kPa  are  close  to  those  used  here.  In  fact, 
Figures  1(e)  and  3(e)  agree  fairly  well  with  Figures  5  and  7  given  by  Edson  and  Roy 
for  their  20  and  200  kPa  experiments  [16]. 
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Figures  1  through  4  provide  information  at  selected  locations  along  the  acoustical 
axis  of  the  array.  These  locations  may  not  provide  the  best  choice  with  respect  to  local 
maxima  in  P+  and  |P_|.  Treating  P+  and  |P_  |  as  a  function  of  a  is  appealing;  however, 
as  Figure  1(b)  suggests,  complicated  constructive  and  destructive  interference  may 
obfuscate  these  waveform  features.  Thus,  a  time-average  energy  flux  density  in  the 
direction  of  propagation  seems  to  provide  a  good  measure.  The  energy  flux  vector 
is  defined  as  j  =  pu,  and  within  the  progressive  wave  approximation  this  reduces  to 
jz  =  p2p2 /poCO)  where  here  p0  is  the  source  amplitude  and  P  is  a  normalized  pressure. 
The  time  average  is  computed  across  the  entire  retarded  time  window, 


0.)  = 


Po 


PqCq{tM  Pn) 


rrM 

Pn)  •'Tm 


P2dr' 


(2) 


such  that  Tm  and  tm  are  the  minimum  and  maximum  of  the  time  window.  Figure 
6  depicts  (jz)  for  the  seven  element  annular  array  (dashed  curve)  and  a  spherically 
focused  piston  transducer  (solid  curve).  As  a  reminder,  the  outer  radius  of  the  piston 
source  and  the  array  are  equivalent,  and  the  relatively  small  separation  between 
adjacent  elements  means  the  piston  transducer  and  array  deposit  the  same  amount 
of  energy  into  the  total  acoustic  field. 


Figure  6.  The  time-averaged  energy  flux  along  the  acoustic  axis  of  the  transducer, 
(a)  po  =  25  kPa;  (b)  p0  =  125  kPa;  (c)  p0  =  225  kPa;  and  (d)  p0  =  325  kPa.  The 
solid  curve  shows  the  results  for  a  spherically  focused  piston  transducer  and  the 
dashed  curve  shows  the  results  for  a  phased  array. 
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Two  important  features  are  immediately  evident  in  Figure  6.  First,  the  curve  for 
the  array  contains  a  strong  peak  near  a  —  0.3.  The  peak  is  due  to  a  combination  of  the 
small  number  of  elements  and  their  finite  dimensions.  The  ratio  of  this  peak’s  value 
to  the  peak  value  at  the  focus  appears  to  be  weakly  dependent  on  po-  This  suggests 
that  a  prefocal  “hot  spot”  may  be  a  feature  of  an  annular  array  source;  its  magnitude 
is  expected  to  diminish  with  an  increasing  number  of  elements.  The  second,  and 
perhaps  most  important,  feature  is  the  significant  difference  between  a  focused  piston 
transducer  and  a  phased  array  at  all  source  amplitudes.  This  difference  is  due  to  the 
finite  dimensions  of  the  individual  array  elements.  With  a  phase  of  ipn  =  GRn,  it  is 
clear  that  the  energy  allotted  to  the  portion  of  the  element  with  R  <  Rn  has  already 
arrived  at  and  passed  through  the  focus,  while  the  energy  allotted  to  the  portion  of 
the  array  element  with  R  >  Rn  has  not  reached  the  focus.  As  the  number  of  array 
elements  is  increased,  the  surface  area  of  each  element  is  reduced,  and  the  phased 
array  result  approaches  the  upper  limit  for  a  spherically  focused  piston  transducer. 
Hence,  any  simulation  of  an  array  of  intense  sources  must  be  carefully  performed. 

It  is  noteworthy  to  briefly  discuss  the  peak  values  of  P+  and  |P_|  although  figures 
for  this  case  study  will  not  be  presented.  At  a  source  amplitude  of  p0  =  25  kPa,  the 
acoustic  fields  from  the  focused  piston  transducer  and  the  array  attain  their  largest 
values  at  a  prefocal  distance.  That  is,  the  geometric  focus  is  at  a  =  1,  but  the  peak 
values  occur  slightly  before  <7  =  1.  As  the  source  amplitude  is  increased,  the  peak 
values  shift  through  the  geometric  focus  and  subsequently  occur  at  a  >  1.  Averkiou 
and  Hamilton  discuss  a  shift  of  peak  values  from  the  geometric  focus  when  P+  and 
|P_|  are  measured  along  the  axis  of  a  focused  circular  piston  operating  at  2.25  MHz 
[17]. 


2.2  Fresh  Water  Lake 

The  computations  discussed  above  in  Section  2.1  and  their  agreement  with  the  re¬ 
sults  of  th  BU  experiments  suggest  that  the  KZK  formalism  is  well  suited  to  analyzing 
nonlinear  wave  propagation  from  a  phased  array.  However,  all  previous  work,  both 
computational  and  experimental,  has  considered  frequencies  and  source  dimensions 
that  are  inappropriate  for  MCM  applications.  This  section  discusses  the  predicted 
nonlinear  behavior  of  a  phased  array  in  a  fresh  water  lake. 

The  characteristics  of  the  simulated  phased  array  have  been  adjusted  to  a  possible 
MCM  system.  The  array  contains  32  elements,  and  the  surface  area  of  each  element 
is  0.0908  m2.  The  innermost  element  has  an  outer  radius  of  0.17  m,  and  the  overall 
outer  radius  of  the  array  is  1.004  m.  The  separation  between  adjacent  elements  is 
0.002  m.  The  source  condition  is  again  given  by  (1),  but  the  frequency  has  been 
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lowered  to  100  kHz  from  the  500  kHz  used  in  the  BU  experiments.  The  sine- wave 
tone  burst  contains  five  cycles.  The  source  amplitude  has  also  been  lowered  to  1,  10, 
20,  and  30  kPa.  A  mine  countermeasure  system  requires  a  standoff  distance  of  tens 
of  meters.  Hence,  the  geometric  focus  of  the  array  has  been  set  to  50  m. 

The  ambient  water  conditions  have  been  adjusted  to  values  that  may  occur  in  a 
fresh  water  littoral  environment.  The  ambient  pressure  is  2  atm,  which  corresponds 
to  an  approximate  depth  of  10  m  (i.e.,  the  array  is  10  m  below  the  lake  surface).  The 
ambient  temperature  and  salinity  are  10°C  and  0  ppt,  respectively.  Hence,  relaxation 
processes  are  neglected  in  this  simulation.  Finally,  the  empirical  rules  for  water  are 
again  used  to  determine  the  parameters  required  for  the  water. 

Figures  7  and  8  compare  retarded  time  traces  at  a  dimensionless  axial  location 
of  a  =  0.8,  where  the  signals  are  near  their  maximum  peak  positive  pressure.  This 
corresponds  to  about  40  m  from  the  array,  and  hence  the  peak  response  occurs  at 
a  prefocal  region  with  respect  to  the  geometric  focus.  Figure  7  demonstrates  that 
linear  superposition  (dashed  curve)  of  the  fields  from  the  array  elements,  when  driven 
independently  but  with  an  appropriate  phase  delay,  fails  to  reproduce  the  acoustic 
field  from  a  properly  phased  array  (solid  curve)  at  high  source  amplitudes.  The  poor 
agreement  shown  in  7(b),  7(c),  and  7(d)  is  due  to  the  neglect  of  the  nonlinear  inter¬ 
action  between  the  fields  from  the  individual  elements.  Figure  8  compares  the  field 
produced  by  the  32  element  annular  array  (dashed  curve)  with  that  produced  by  a 
focused  piston  transducer  (solid  curve)  with  comparable  source  parameters.  Within 
the  resolution  of  Figure  8,  these  results  are  essentially  indistinguishable.  Hence,  the 
claim  in  Section  2.1  that  better  agreement  between  an  annular  array  and  focused  pis¬ 
ton  transducer  will  occur  as  the  number  of  array  elements  increases  is  substantiated. 
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Figure  7.  Linear  superposition  of  the  nonlinear  fields  from  the  independent  array- 
elements  (dashed  curve)  and  the  acoustic  field  from  a  properly  phased  array  (solid 
curve),  (a)  p0  =  1  kPa,  (b)  p0  =  10  kPa,  (c)  p0  =  20  kPa,  and  (d)  p0  =  30  kPa. 
The  dimensionless  axial  location  is  a  =  0.8. 


Figure  8.  The  acoustic  fields  from  a  properly  phased  array  with  32  elements 
(dashed  curve)  and  a  focused  piston  source  (solid  curve).  The  sources  produce 
similar  results.  The  source  amplitudes  are  (a)  1,  (b)  10,  (c)  20,  and  (d)  30  kPa.  The 
dimensionless  axial  location  is  a  =  0.8. 
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A  comparison  of  the  time-averaged  energy  flux  density,  (jz),  is  depicted  in  Figure  9. 
The  dashed  curve  is  the  result  obtained  for  the  array,  and  the  solid  curve  is  the  result 
produced  by  the  focused  piston  source.  Two  features  are  immediately  evident.  First, 
the  peak  occurs  near  cr  =  0.7  for  p0  =  1,  10,  and  20  kPa,  with  a  slight  inward  shift  for 
p0  =  30  kPa.  The  prefocal  maximum  is  similar  to  results  reported  by  Hutchins  et  al. 
[18]  for  annular  arrays  and  spherical  bowl  sources  under  linear  acoustic  conditions. 
However,  the  additional  shift  for  p0  =  30  kPa  is  due  to  the  nonlinearity  of  the  fluid. 
A  second  feature  is  the  narrowing  of  the  main  peak  at  the  higher  source  amplitudes 
of  20  and  30  kPa.  This  suggests  the  focal  volume  is  slightly  compressed  at  the  higher 
source  amplitudes,  and  the  compression  localizes  the  energy  to  a  smaller  volume 
about  the  focus.  To  fully  quantify  the  acoustic  field  in  the  vicinity  of  the  focus,  (jz) 
should  be  computed  at  off-axis  locations  along  an  axis  transverse  to  the  ^-direction. 
This  procedure  has  not  been  implemented;  however,  the  peak  positive  pressure  as  a 
function  of  the  radial  coordinate  at  the  focus  (i.e.,  the  beam  width)  has  been  found 
to  be  in  agreement  with  linear  acoustics.  Finally,  under  ordinary  linear  acoustics 
when  the  source  amplitude  is  doubled  (or  tripled),  then  the  resulting  change  in  (jz)  is 
essentially  a  factor  of  4  (or  9).  Inspection  of  Figure  9  shows  that  a  similar  result  does 
not  hold.  This  behavior  is  due  to  the  repartitioning  of  energy  from  the  fundamental 
frequency  to  its  higher  harmonics  and  hence  the  subsequent  higher  energy  loss  from 
thermoviscous  absorption. 


Figure  9.  Comparison  of  the  time-averaged  energy  flux  density  for  a  properly 
phased  array  with  32  elements  (dashed  curve)  and  a  focused  piston  transducer  (solid 
curve).  The  results  agree,  (a)  po  =  1  kPa,  (b)  po  =  10  kPa,  (c)  po  —  20  kPa,  and 
(d)  po  =  30  kPa.  The  dimensionless  axial  location  is  a  =  0.8. 
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3.  THEORETICAL  DEVELOPMENT 


3.1  Augmented  Khokhlov-Zabolotskaya-Kuznetsov  Equation 


Cleveland  et  al.  discuss  a  KZK  nonlinear  parabolic  wave  equation  which  has  been 
augmented  by  including  the  relaxation  processes  within  the  fluid  [13].  One  form  for 
the  augmented  KZK  equation  for  pressure  fluctuation,  p,  is 


d2p  _  Cp  2  6  &p  (3  d2p2  ^  1  .  d3p 

dt'dz  2  ±P  2 cldt'3  2p0Co  dt'2  „  2 pQcl  dt'3  ’ 


where  ±  denotes  directions  transverse  to  the  propagation  direction  z  [19],  f  =  t  —  z/co 
is  a  retarded  time  (that  facilitates  defining  a  temporal  window  which  moves  with  a 
pulse),  t  is  time,  and  <5  is  the  diffusivity  of  sound;  <5  =  [C+4r?/3+K(c“1  —  c~ 1)]/po-  The 
bulk  and  shear  viscosity  of  the  fluid  are  (  and  rj,  k  is  the  thermal  conductivity,  c*,  is 
the  specific  heat  at  constant  volume,  and  Cp  is  the  specific  heat  at  constant  pressure. 
The  relaxation  operator  for  the  uth  relaxation  process  is 


Ipv  =  PoC2o 


mvtv 
1  4*  tv  dt1 


d  ’ 


(4) 


tv  is  a  relaxation  time,  and  mu  is  a  dimensionless  quantity  characterizing  the  disper¬ 
sion  of  the  relaxation  process  [20]. 


Although  a  progressive  wave  assumption,  p  =  PqCqUz ,  is  required  to  arrive  at  (3), 
the  first  term  on  the  right-hand  side  accounts  for  diffraction.  That  is,  some  energy 
propagates  in  the  transverse  direction,  whereas  it  has  been  assumed  that  propagation 
primarily  occurs  along  the  z  axis.  The  second  term  accounts  for  thermoviscous  ab¬ 
sorption.  It  does  not  include  relaxation  processes,  which  are  explicitly  included  by  the 
last  term.  The  third  term  represents  quadratic  nonlinearity.  Harmonic  generation, 
leading  to  shock  formation,  is  accounted  for  by  this  term.  Finally,  the  last  term  in 
(3)  is  a  summation  over  all  relaxation  processes  in  the  fluid.  For  clean  fresh  water, 
i/v  =  0,  and  (3)  reduces  to  the  standard  KZK  equation  [21]. 


3.2  Axisymmetric  Circular  Source 

When  the  source  condition  in  the  z  =  0  plane  corresponds  to  an  axisymmetric, 
circular  region  such  as  a  piston  transducer  or  annular  array,  then  (3)  can  be  written 
as 

a2p  _  co  (cP_  ld_\  _5_&p  /3  02p2 

dt'dz  2  \dr 2  r  dr )  P  +  2cq  dt 13  2poCo  dt'2 
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(5) 


where  r  is  the  radial  coordinate.  There  are  no  known  analytic  solutions  to  (5),  so  it 
must  be  solved  numerically.  It  is  convenient  to  reduce  this  equation  to  a  dimensionless 
form.  Introduce  the  following  scaling  transformations: 


P'=p/Po,  T'  =  u0t',  R'  =  r/a, 
d  d 
dt^^d?’ 

d2  Id  1  /  d2  1 

dr2  +  ~r~dr  a2  \&R'2  +  R'  dR' )  ’ 


(6) 

(7) 

(8) 


where  w0  is  a  characteristic  angular  frequency  of  the  source  waveform,  and  a  is  the 
radius  of  the  circular  source.  These  transformations  lead  to 

d2P'  _  1  (  d2  ±_d_\  ,  g3p/  !  g2 pg 

dr'dz  ~  4zT  \dR'2  +  R' dR' )  a°  dr'3  2zs  dr'2 

muljJo  (  Tv  \  dsP 

+  ?  IcT  vT+^J- j  Hr*'  (9) 

Equation  (9)  introduces  three  new  constants  that  characterize  an  axisymmetric  cir¬ 
cular  source  and  the  nonlinear  wave  propagation.  The  first,  zr  =  w0a2/2c0,  is  the 
well-known  Rayleigh  distance  for  a  circular  aperture.  This  distance  is  approximately 
the  axial  location  where  a  linear  acoustic  field  transitions  from  the  so-called  nearfield 
to  the  farfield.  The  second,  o0  =  ^o/2co»  is  the  thermoviscous  attenuation  coeffi¬ 
cient.  For  an  inhomogeneous  plane  wave  and  linear  acoustics,  ao  is  the  distance  the 
plane  wave  must  travel  for  it  to  attenuate  to  e_1  of  its  original  amplitude.  The  third 
constant,  =  p0cl//3uj0po,  is  the  plane-wave  shock  formation  distance.  This  is  the 
distance  an  initially  sinusoidal  disturbance  must  propagate  in  order  for  the  cumula¬ 
tive  effects  of  the  quadratic  nonlinearity  to  produce  shock  formation.  In  addition,  the 
dimensionless  relaxation  time  is  tv  =  uJotu. 


All  quantities  in  (9)  have  been  scaled  by  an  appropriate  value  except  the  z  coor¬ 
dinate.  Previous  researchers  have  used  either  the  Rayleigh  distance  for  a  flat  piston 
source  or  the  focal  length  of  a  weakly  focused  source  for  this  scaling.  A  third  possible 
choice  for  the  length  scale  is  the  shock  formation  distance.  To  accommodate  each  of 
these  choices,  an  arbitrary  length  scale  l  will  be  introduced  such  that 


t _ z  d_  1  d_ 

°  V  dz~*  l  da'  ’ 


(10) 


and  (9)  becomes 

d2P'  l  (  d2  1  d  \  ,  d3P'  l  d2P'2 

dr' da'  ~  4 zT  \dR'2  +  R' dR' )  a°  dr'3  2zs  dr'2 
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+ 


771i/OJqI  (  Tv 
»  2Co  Vl+T,,^ 


dr'3  ’ 


(11) 


For  a  given  pulse  excitation  in  the  a'  =  0  plane  as  a  function  of  r'  and  R',  (11)  can 
be  solved  numerically  by  finite-difference  and  operator-splitting  methods  (see  below). 


Previous  researchers  have  found  that  an  additional  coordinate  transformation  is 
useful,  especially  when  diffraction  becomes  important.  This  change  of  coordinates  is 
sometimes  referred  to  as  a  “stretched  coordinate”  transformation  because  a  rectan¬ 
gular  grid  in  the  transformed  coordinate  space  (a,  R )  corresponds  to  a  grid  that  has 
been  stretched  in  the  original  coordinate  space  (a1,  R!).  A  complete  discussion  of  this 
transformation  is  given  by  Lee  [22].  A  convenient  form  for  the  “stretched  coordinate” 
transformation  is 


R! 


a  —  a 


P=  (!  +  {✓)*>', 


T  — 


£R'2 
1  +  £<r' 


(12) 


where  the  significance  of  the  parameter  £  s  discussed  below.  The  differential  operators 
in  (11)  become 


9 2 

dR'2  + 


d_  d_ 

dr'  dr  ’ 

A  _»  a  eg  g  |  — 

da'  da  1  +  ZadR  dr ’ 

Id  1  (  d2  1  d  \ 

R'  dR'  (1  +  £<r)2  [dR2  +  RdR) 

1  +  tadRdr^  ?  dr2 


4£  d 
1  +£a  dr 


(13) 

(14) 


(15) 


Substitution  of  (12)  to  (15)  into  (11),  integration  with  respect  to  r',  and  some  tedious 
algebraic  manipulations  give 


dP 

da 


l 


+ 


+ 


4*r(l  +  &)2 
l 


L( 


8 2  1  9\„;,  ,d2P 

gm  +  RdR  )pdT  +aolarn 


_ pdP  sp  m^oJol 

zs(l  +  £cr)  dr  V  2co 
£  A  l\  d(RP) 


1  +  T„ 


dr  - 


1  +£<7 


H) 


2  d2 


dR 


Z  R 


'i-T 


d2p 
dr2 
\  dP 
I  dr' 


(16) 


If  £  =  0,  then  (16)  reduces  to  essentially  an  identity  transformation  where  the  primes 
in  (11)  are  dropped.  When  diffraction  effects  are  expected  to  be  minimal,  as  with 
a  weakly  focused  source,  then  (11)  can  be  used  directly  and  l  is  the  focal  length.  If 
£  =  i  and  l  =  zT,  the  transformation  leads  to  the  “stretched  coordinate”  system. 
Equation  (16)  is  appropriate  when  diffraction  is  expected  to  be  important.  Finally, 
the  last  two  terms  in  (16)  can  be  ignored  under  these  conditions  on  £  and  l. 
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3.3  Rectangular  Source 


When  the  source  is  rectangular  with  dimensions  x0  and  y0,  (3)  can  be  written  in 
the  following  dimensionless  form: 


dP_ 

da 


+ 


Axr 


d2P 


l 


»  +  4ft 


f 

J-oo 


d2PJ  ,  ,d2P  ,  l  U8P 

—-dr  +  a^l-T - P— 

dY2  dr 2  zs  dT 


cPP_ 

dT2 


(17) 


The  distances  xr  =  cu0a^/2co  and  yr  =  w0^/2c0  are  defined  in  a  manner  consis¬ 
tent  with  the  definition  of  the  Rayleigh  distance  zr,  and  X  and  Y  are  dimensionless 
Cartesian  coordinates.  Although  the  terms  accounting  for  thermoviscous  absorption, 
quadratic  nonlinearity,  and  relaxation  processes  in  (17)  are  similar  to  those  terms 
in  (16),  these  terms  have  implicit  dependence  on  the  Cartesian  coordinates  through 
P(X,  Y,  a,  t).  A  transformation  to  a  “stretch  coordinate  system”  could  be  performed, 
but  it  has  not  been  pursued.  Baker  et  al.  describe  a  transformed  beam  equation  for 
a  rectangular  aperture  that  is  based  on  a  standard  KZK  equation  [23].  Their  result 
is  essentially  the  “stretched  coordinate”  transform  of  (17)  without  relaxation. 
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4.  NUMERICAL  SOLUTION 


There  are  no  known  analytic  solutions  to  (3).  Thus  numerical  solutions  are  the 
only  available  means  to  investigate  nonlinear  wave  propagation  from  a  bounded  beam 
source.  For  a  given  pulse  source  condition  F(R,a  =  0  ,r),  the  axisymmetric  source 
in  Section  3.2  is  a  two-dimensional  problem  in  ( R ,  r)  of  forward  marching  from  a  to 
<7  +  Aa.  If  A  a  is  small,  then  (16)  can  be  solved  via  an  operator-splitting  method  [24]. 
That  is,  the  following  system  of  equations  is  equivalent  to  (16): 


dP  1  fT  (  92  \  1  9\  Pd-' 

da  ~  4G(l  +  Za)2J-oo\dR?  RdRj  T 

dP_  _  Ad2P 

da  ~  Adr 2’ 
dP  _  NP  dP 

da  1  +  £a  dr  ’ 


(18) 

(19) 

(20) 

(21) 


(subject  to  the  above  conditions  on  £  and  /).  In  these  equations,  G  =  zr/l  is  known 
as  the  focusing  gain,  A  —  a0l  is  a  thermoviscous  absorption  constant,  N  =  l/zs  is  a 
nonlinearity  constant,  and  Cv  =  mJ/rI/w0i/2co  is  a  dimensionless  parameter.  Equation 
(21)  is  applied  to  each  relaxation  process. 


Numerical  integration  of  (18)— (21)  requires  the  specification  of  a  discrete,  finite 
domain  in  (R,  r)  and  appropriate  conditions  at  the  edge  of  this  domain.  The  finite 
domain  is  R  £  [0,  RM]  and  r  £  [rm,rM];  RM  is  chosen  large  enough  to  eliminate  a 
possible  reflection  from  the  nonabsorbing  boundary.  The  minimum  and  maximum 
dimensionless  retarded  time,  rm  and  tm,  are  chosen  such  that  a  zero-amplitude  pe¬ 
riod  of  time  precedes  and  follows  the  finite  pulse,  permitting  the  pulse  to  propagate 
without  edge  effects  from  the  finite  domain.  Furthermore,  the  following  conditions 
are  imposed: 


dP 

P(a,R,Tm)  =  0,  P(a,  R,tm)  =  0,  P{a,  Rm,t)  =  0,  — 


0. 


(22) 


R= 0 


The  last  condition  is  essentially  a  symmetry  condition  along  the  axis  of  the  circular 
source.  Introduction  of  the  fixed  increments  At  and  A R  allows  a  discrete  grid  to  be 
defined  such  that  r  £  [ Tm,TM ]  <=>  i  €  [0,/]  and  R  £  [0,  Rm]  <*=»  j  £  [0,  J],  where 
P(a,  R,  t)  =  P(ak,jAR,  Tm+?Ar)  =  PtA  and  the  first  three  conditions  in  (22)  reduce 
to  P*j  =  0,  Pfd  =  0,  and  P£j  =  0. 


Standard  finite-difference  (FD)  techniques  can  be  implemented  to  solve  (18),  (19), 
and  (21).  The  integration  in  (18)  is  first  approximated  with  a  standard  trapezoid 
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quadrature.  The  details  of  reducing  these  equations  to  FD  algorithms  are  given  by 
Lee  [22]  and  therefore  are  omitted  from  the  present  report.  Under  source  conditions 
of  interest  to  Lee,  it  was  observed  that  implicit  backward  FD  algorithms  were  re¬ 
quired  near  the  source  plane.  At  a  sufficient  distance  from  this  plane,  the  numerical 
integrations  could  be  switched  to  Crank-Nicolson  algorithms.  The  numerical  advan¬ 
tages  of  the  Crank-Nicolson  scheme  with  respect  to  an  implicit  backward  method  are 
discussed  at  length  in  standard  texts  on  FD  techniques  [25].  Finally,  analyses  of  the 
truncation  errors  show  that  the  following  algorithms  are  accurate  in  A  a ,  AR,  and 
At  to  at  least  first  order. 

The  equation  governing  nonlinear  steepening,  (20),  yields  an  analytic  solution.  If 
the  dimensionless  pressure  is  P(a,R,r),  then  it  can  be  shown  by  direct  substitution 
that 

P(a  +  Aa,R,r)  =  P(a,R,r  +  NP<j>a)  (23) 

is  a  solution  to  (20).  The  factor  4>a  in  our  notation  is 


Again,  the  details  are  given  by  Lee  [22],  The  numerical  algorithm  is  described  below. 
Furthermore,  a  single  algorithm  suffices  through  the  entire  a  domain. 


4.1  Diffraction  Algorithm 


The  diffraction  integral,  the  symmetry  condition  along  the  axis  of  the  source, 
and  the  nonabsorbing  boundary  condition  at  R  —  Rm  suggest  that  the  FD  algorithm 
should  give  explicit  expressions  for  j  =  0, 1  <  j  <  J— 1,  and  j  =  J— 1  for  i  €  [1,7—1]. 
The  implicit  backward  FD  (IBFD)  algorithm  is 


(1  -  D2)P$l  +  D2Pff 
Dsf-P, +  (1  -  Dt)Pft'  + 

D,fJ-iP?jU  +  (1  -  Di)P&h 


Pl 0  -  A(Si  -  So), 

(25) 

P?4  +  AS, 

D,{f+Sj+ 1  +  /-S,_,), 

(26) 

Pl,  ,  -  A/7_,Sj-2 

D2SJ-I, 

(27) 

where 


Dn  =  - 


ArAjfeCr 


i— 1 


nG(AR)2(l  +  £<7fc+i)2  ’ 


mj  * 


m—  1 


(28) 
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The  first  equation  of  (28)  is  a  constant  for  a  given  A *<?  step.  The  subscript  k  on 
A  indicates  that  an  adaptive  step-size  is  possible  in  a.  In  fact,  in  the  vicinity  of  a 
strong,  rapid,  change  in  pressure,  Acr  needs  to  be  small  to  prevent  the  occurrence 
of  an  unphysical  multivalued  waveform  during  the  nonlinear  step  described  by  (20). 
The  second  and  third  equations  in  (28)  arise  from  the  diffraction  integral  in  (18)  and 
its  subsequent  approximation  by  trapezoid  quadrature.  (Note:  The  values  ■  in  the 
summation  occur  at  previous  moments  in  time.)  The  Crank-Nicolson  FD  (CNFD) 
algorithm  is 

i(l  -  DJQft1  +  D,Q*t'  =  1%  -  Dt(Si  -  S'„),  (29) 

D32KQ%l1  +  l(l-D,)Q'Ql  +  D32f+Ql8ll  =  Pfj  +  DsS' 

-  DM+S'j+l  +  (30) 

Dnfi-lQ&U  +  ^(1  -  AOQtJi,  =  Kl- 1  ~  Ae/7-iSJ-a 

+  A>SJ-i,  (31) 

where 

Qt y  =  Pi/1  +  Plj.  ^  =  E  Q£j-  (32) 

m=  1 

Inspection  of  (25)-(27)  and  (29)-(31)  shows  that  the  IBFD  and  CNFD  algorithms 
form  tridiagonal  systems  of  equations.  A  standard  matrix  method,  known  as  the 
Thomas  algorithm,  permits  a  rapid  solution  to  these  equations. 


4.2  Thermoviscous  Absorption  Algorithm 


The  IBFD  algorithm  for  the  equation  governing  thermoviscous  absorption  can  be 
written  as 

A'Pffj  +  (1  -  2A')P?f  +  A'P&'j  =  (33) 

where  A!  =  —  AAjtcr/(Ar)2.  The  corresponding  CNFD  algorithm  is 


Y Pffj  +  (1  -  A ')/?/'  +  Y Pi/j  =  (1  +  A')P?J  -  Y (Pt-u  +  Pwj)-  (34) 


When  edge  conditions  Pfij  and  Pjj  are  introduced  into  (33)  and  (34),  it  is  seen  that 
tridiagonal  systems  of  equations  are  recovered. 
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4.3  Quadratic  Nonlinearity  Algorithm 


As  stated  previously,  the  quadratic  nonlinearity  manifests  an  analytic  solution. 
The  numerical  implementation  of  (23)  and  (24)  involves  two  steps.  First,  the  wave¬ 
form  for  each  radial  value  j  is  distorted  onto  a  nonuniform  f  grid.  This  is  accomplished 
by 


=  Tti  ~  NPtjc 


(35) 


Here,  rtk  are  the  uniformly  spaced  dimensionless  retarded  time  points  such  that  At  = 
rfc.  _  and  f fj  are  the  new  retarded  time  points  which  are  not  uniformly  spaced. 
Both  Pf-  and  ^’are  evaluated  at  the  current  ak  location.  As  can  be  seen  in  Figure 
10,  this  phase  distortion  “slides”  the  pressure  amplitude  values  to  new  retarded  times. 


is  the  final  value  of  pressure,  Pkfl,  after  the  resampling  step. 


The  second  step  entails  a  resampling  of  the  distorted  waveform  back  onto  the  uniform 
retarded  time  grid.  This  step  is  necessary  because  the  finite-difference  algorithms  for 
diffraction,  thermoviscous  absorption,  and  relaxation  assume  a  uniform  grid.  Simple 
linear  interpolation  is  used  to  resample  the  waveform; 


pfc+1  __ 
hj 


pk 

rm+l,j 


Tk 

/m+l  J 


(t*,  -  fij) + p, 


k 


(36) 
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where  is  the  pressure  amplitude  corresponding  to  f£j.  A  pictorial  description 
of  these  steps  is  contained  in  Figure  10. 

Physically,  a  valid  solution  of  (23)  is  constrained  to  be  a  single-valued  function  of 
retarded  time.  The  algorithm  here  does  not  guarantee  this  requirement,  and  hence 
application  of  (35)  means  that  a  strict  ordering  of  <  f£+lj  must  be  maintained. 
The  desired  ordering  can  be  imposed  by  selecting  to  be  an  appropriately  small 
value;  in  fact,  an  adaptive  step-size  in  o  can  be  implemented  by  requiring 


Afccr  < 


(37) 


where  the  subscript  “max”  denotes  the  maximum  slope  in  the  newly  distorted  wave¬ 
form. 


As  a  final  note,  the  nonlinear  steepening  of  a  waveform  shifts  energy  from  lower- 
frequency  components  to  higher  frequencies.  The  algorithm  in  (35)  is  essentially  exact 
(to  within  numerical  round-off  errors)  and  is  certainly  more  accurate  than  the  finite- 
difference  algorithms.  However,  the  resampling  step  produces  an  artificial  smoothing 
of  the  waveform  because  At  is  always  larger  than  some  step  sizes  on  the  nonuniform 
grid.  The  smoothing  acts  as  a  low-pass  filter.  Hence  for  accurate  numerical 
simulations  At  must  be  chosen  to  accommodate  the  highest  frequency  of  interest. 


4.4  Relaxation  Algorithm 

The  IBFD  algorithm  for  the  ut h  relaxation  process  is 

-KP&J  +  (1  +  2 Cljpft'  -  KP&'j  =  Pi,i  +  <(P?+1J  -  Pt-u).  (38) 

where  C'v  —  O,  (Act)* /(At)2,  t'  =  t„/2At,  and  A  *  =  C"  ±  t'.  The  Crank-Nicolson 
method  yields 

-A tP?-\]j  +  2(1  +  Cl)P?f  -  A ;Pg?j  =  2(1  -  CJPfj  +  KPUi  +  KPLj,  (39) 

such  that  A/  =  ±  t'. 
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APPENDIX 

Empirical  Material  Parameters  for  Water 

In  Section  3,  several  symbols  are  introduced  to  represent  physical  properties  of 
the  fluid.  Standard  handbooks  provide  limited  tables  for  these  properties  for  water; 
however  many  are  restricted  to  small  ranges  of  temperature,  ambient  pressure,  and 
salinity  [26].  One  is  forced  to  scour  the  literature  to  obtain  a  complete  set  of  param¬ 
eter  for  any  given  computation.  The  purpose  of  this  appendix  is  to  compile  several 
empirical  relationships  for  these  parameters.  An  empirical  rule  for  the  parameter  of 
nonlinearity,  B /A,  for  water  is  not  available,  so  a  table  of  values  from  the  literature 
is  provided.  The  table  permits  linear  interpolation  for  values  not  contained  explicitly 
in  the  table. 

The  small-signal  sound  speed,  c0,  can  be  approximated  from  the  Chen-Millero- 
Li  equation  [27],  which  is  an  empirical  fit  to  a  large  range  of  data.  It  has  explicit 
dependence  on  temperature,  ambient  pressure  relative  to  1  atm,  and  salinity.  These 
quantities  are  denoted  respectively  by  T,  P,  and  S.  The  Chen-Millero-Li  equation  is 

Co  =  cw(P ,  T)  -  cc(P,  T )  +  a(P,  T)S  +  b(P,  T)S3/ 2  +  d(P)S2,  (40) 

where  the  auxiliary  functions  are 

cw(P,T )  =  1402.388  +  5.03711T  -  0.0580852T2  +  3.3420 xlO"4T3 

_  1.4780  xlO"6T4  +  3.1464  xHT9T5  +  (0.153563  +  6.8982  xlO"4T 

-  8.1788  x  10~6T2  +  1.3621  x  10-7r3  -  6.1185  x  10"10T4)P 

+  (3.1260  xlO"5  -  1.7107 x  10-6T  +  2.5974  x  10-8T2  -  2.5335  x  10-10T3 

+  1.0405  xlO"12T4)P2  -  (9.7729  xlO-9  -  3.8504 x  10"10T 

+  2.3643  x  10_12T2)P3,  (41) 

cc(P,T)  =  (0.0029  -  2.19 xl0_4T+ 1.4 xlO_5T2)P- (4.76 xlO-6 

-  3.47 xl0_7r  + 2.59 x10_8T2)P2  + 2.68  xlO_9P3,  (42) 

a(P,T)  =  1.389  —  0.01262T  +  7.164  x  10_5T2  +  2.006x  10_6T3  —  3.21  x  10_8T4 

+  (9.4742  xlO"5  -  1.2580  xl0-5T- 6.4885  xlO-8P2  +  1.0507  xlO-8P3 

-  2.0122  x  10-loT4)P  -  (3.9064  x  10"7  -  9.1041  x  10“9T 

+  1.6002  x  10_1°T2  —  7.988  xl0_12T3)P2  +  (l.lOOxlO-10 

+  6.649  x  10_12T  —  3.389  x  10_13T2)P3,  (43) 

b(P,T)  =  —0.01922  —  4.42  x  10-5T  +  (7.3637  x  10-5  +  1.7945  x  10~7T)P ,  (44) 

d(P)  =  0.001727-  7.9836  xlO~6P.  (45) 

Equations  (40)  to  (45)  are  subject  to  the  following  restrictions:  0  <  T  <  40°C, 
0  <  S  <  40  ppt,  and  0  <  P  <  1000  bar.  Within  the  applicable  ranges  on  T,  P,  and 
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S,  the  Chen-Millero-Li  formula  is  claimed  to  be  accurate  to  ±0.05  m/s,  where  Cq  has 
units  of  meters  per  second.  Finally,  it  is  restated  that  P  has  units  of  bar. 

The  density  of  the  ambient  fluid,  specific  heat  at  constant  pressure,  specific  heat 
at  constant  volume,  and  the  thermal  conductivity  are  determined  by  the  following 
empirical  equations  [28]: 


Po 

Cp 

Cy 

K 


j  999.7  +  0.048  x  10"5P  -  0.088AT  -  0.007(AT)2, 

\  1027  +  0.043  x  10~5P  -  0.16AT  -  0.004(AT)2  +  0.75AP, 

j  4192  —  0.40  x  10~5P  —  1.6AT, 

\  3988  -  0.23  x  10"5P  +  0.54AT  -  5.4AS, 

J  Cp[l  -  0.0011(1  +  AT/6  +  0.0024  x  10"5P)2], 

\  cp[ 1  -  0.0041(1  +  AT/20  +  0.0012  x  10'5P  +  0.012AS)2], 

0.597  +  0.0017AT  -  7.5  x  10~6(  AT)2. 


(46) 

(47) 

(48) 

(49) 


The  upper  expressions  for  p0,  Cp,  and  cv  apply  to  fresh  water,  and  the  lower  expressions 
apply  to  sea  water.  The  salinity  is  taken  relative  to  35  ppt  (i.e,  A S  =  S  —  35),  and 
the  relative  temperature  is  AT  =  T-  273.16  K.  Unlike  the  Chen-Millero-Li  equation, 
P  is  absolute  pressure  in  units  of  pascal  and  T  is  specified  in  kelvin.  The  units  are 
kilograms  per  cubic  meter  for  po,  joules  per  kilogram  per  kelvin  for  Cp  and  cv,  and 
watts  per  meter  per  kelvin  for  k.  Finally,  Pierce  states  that  the  thermal  conductivity 
is  essentially  constant  with  respect  to  changes  in  salinity  and  pressure.  Hence,  (49)  is 
sufficient  for  both  fresh  and  sea  water,  and  it  should  deviate  by  no  more  than  a  few 
percent  from  the  actual  value  for  k. 

The  shear  viscosity  of  water  appears  to  be  insensitive  to  changes  in  salinity  and 
ambient  pressure  but  varies  with  temperature.  A  standard  handbook  [26]  containing 
physical  properties  for  water  gives  the  following  empirical  formulae: 


log  77  = 


_ 1301 _ o  Qf)933  n  <  T  <  2D 

998. 333+8. 1855AT+0.00585(  AT)2  '3'3U/'3,3’  U  ±  1  ±  ^u> 

_  1.3272AT+aO01053(AT)2  +  g  6772  x  1Q-3  20  <  T  <  100, 

Ai  +125 


where  temperature  is  given  in  Celsius  such  that  AT  =  T  -  20.  The  bulk  viscosity  of  a 
fluid  is  often  ignored  because  of  assumptions  of  irrotational  flow  and  incompressible 
fluid  dynamics.  Pierce  states  that  ( is  relatively  constant  for  changes  in  P  and  S  and  it 
has  a  weak  dependence  on  temperature.  Bulk  viscosity  for  water  can  be  approximated 
by 

C/77  = -0.29T/60.0  + 3.01,  (51) 

where  the  temperature  range  is  0  <  T  <  60°C.  The  units  for  77  and  C  are  centipoise. 
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The  parameter  of  nonlinearity,  B/A ,  is  related  to  the  coefficient  of  nonlinearity 
by  -  1  +  B /2A  for  a  fluid.  Limited  data  inhibit  an  accurate  empirical  formula 
to  estimate  B/A  for  either  fresh  water  or  sea  water  [29,  30].  Estimates  of  B/A  are 
obtained  via  linear  interpolation  from  the  values  listed  in  Table  1.  Pierce  states 
that  nearly  99.5%  of  sea  water  falls  into  a  salinity  range  from  33  to  37  ppt.  Hence, 
estimating  B/A  from  this  table  is  adequate  for  the  numerical  computations. 


Table  1.  B/A  values  for  fresh  and  sea  water. 


T  (°C) 

0  10  20  30  40  60  80  100 

5  —  0 

5  =  33 
5  =  35 

4.2  5.0  5.4  5.7  6.1  6.1 

4.9  5.1  5.2  5.4 

5.25 

The  empirical  models  for  sea  water  are  completed  once  formulae  for  the  di¬ 
mensionless  relaxation  time,  tu  =  LOoty,  and  dimensionless  dispersion  parameter, 
Cv  =  ml/riyujol/2co,  are  determined  for  the  various  relaxation  processes.  As  stated 
in  Section  1,  the  presence  of  H3B03  and  MgS04  in  sea  water  provides  the  dominate 
energy  loss  through  relaxation  mechanisms.  Francois  and  Garrison  give  the  following 
empirical  formula  for  the  attenuation  coefficients: 


a 


V 


AvPvfvp 

P  +  P 


(52) 


where  fu  is  the  relaxation  frequency  and  Av  and  Pv  are  coefficients  that  depend  on 
P,  S,  and  T  as  well  as  the  pH  for  H3B03  [1,  2].  Hence  tv  and  mu  can  be  related  to 
av  via  consideration  of  (21)  and  an  inhomogeneous  plane- wave  exp {ikz  -  iut),  where 
au  =  ^(k).  Simple  algebraic  manipulations  lead  to 

tv  —  l/2nfv,  ttiu  —  AVP (b3) 


Ty  =  Uj/2'Kfv,  Cy  =  (1  +  Ty)(Xyl.  (54) 

The  units  of  fu  and  Av  are,  respectively,  hertz  and  nepers  per  meter  per  hertz  while 
Py  is  unitless.  For  H3B03,  the  parameters  are 

/x  =  2800(5/35)1/210(4_1245/5),  (55) 

Ai  =  1.0200  x  10(0  78pi?_n)/c,  (56) 

Px  =  1,  (57) 

c  =  1412  +  3.21T  +  1.195  +  0.0I67D,  (58) 
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where  0  =  T  +  273  and  T,  S,  and  D  (depth)  have  units  of  Celsius,  parts  per  thousand, 
and  meters,  respectively.  For  MgS04,  these  parameters  are  identified  as 


f2  =  [8.17xlO(11-1990/e)]/[l  +  0.0018(5  -35)], 
A2  =  2.4683  x  10_65(1  +  0.0025T)/c, 

P2  =  l-(1.37xl0-4- 6.2x10 ~9D)D, 

where  c  in  (60)  is  given  by  (58). 


(59) 

(60) 
(61) 
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